Single level regression

Considérons le jeu de données suivant, contenant une variable catégorielle et deux variables quantitatives, et concernant le QI d’enfants et celui de leur mère ainsi que le niveau de formation de celles-ci, supérieur ou non.

kidiq <- read_dta('data/child.iq/kidiq.dta')
head(kidiq)

## # A tibble: 6 × 5
##   kid_score mom_hs mom_iq mom_work mom_age
##       <dbl>  <dbl>  <dbl>    <dbl>   <dbl>
## 1        65      1  121.         4      27
## 2        98      1   89.4        4      25
## 3        85      1  115.         4      27
## 4        83      1   99.4        3      25
## 5       115      1   92.7        4      27
## 6        98      0  108.         1      18

Modèles sans interactions

Nous voulons expliquer le QI des enfants au moyen des deux variables concernant leurs mères. Voici un premier modèle explicatif :

attach(kidiq)
fit.3 <- lm (kid_score ~ mom_hs + mom_iq)
display(fit.3)

## lm(formula = kid_score ~ mom_hs + mom_iq)
##             coef.est coef.se
## (Intercept) 25.73     5.88  
## mom_hs       5.95     2.21  
## mom_iq       0.56     0.06  
## ---
## n = 434, k = 3
## residual sd = 18.14, R-Squared = 0.21

—- Collapsible section —-

Explications : les coefficients et leur variance

Si \(H\) est la matrice de projection sur \(\text{Esp}(X)\) l’espace engendré par les vecteurs colonnes de \(X\), alors

\[ Y = H\cdot Y+(I-H)\cdot Y \]

et tout vecteur \(v\) de cet espace s’écrit comme une combinaison linéaire \(v=\alpha\cdot X\), si bien que

\[ \begin{aligned} \left<X\cdot\alpha\,|\, (I-H)\cdot Y\right> &= 0 \quad, \forall\alpha\\ \alpha^T\cdot X^T\cdot(I-H)\cdot Y &= 0 \quad\text{(simplifiable par }\alpha \text{ puisque }\forall\alpha) \\ X^T\cdot Y &= X^T\cdot H\cdot Y \\ \text{avec}\\ H\cdot Y &= X\cdot\hat{\beta} \\ \text{si bien que} \\ \hat{\beta} &= (X^T\cdot X)^{-1}\cdot X^T\cdot Y \\ \text{et} \\ H &= X\cdot(X^T\cdot X)^{-1}\cdot X^T \\ \end{aligned} \]

\(\hat{Y}=H\cdot Y\) est la projection de \(Y\) sur \(\text{Esp}(X)\).

\[ \begin{aligned} Var(\hat{\beta}) &= Var((X^T\cdot X)^{-1}\cdot X^T\cdot Y) \\ &= (X^T\cdot X)^{-1}\cdot X^T\cdot \underbrace{Var(Y)}_{=Var(\varepsilon)\\=\sigma^2\cdot I}\cdot \left[(X^T\cdot X)^{-1}\cdot X^T\right]^T \\ (\text{puisqu'on a en général que } & Var(M\cdot Z)=M\cdot Var(Z)\cdot M^T)\\ &= \sigma^2\cdot(X^T\cdot X)^{-1} \end{aligned} \]

On obtient un estimateur de \(Var(\hat{\beta})\) en replaçant dans l’expression obtenue ci-dessus \(\sigma^2\) par son estimateur sans biais \(\hat{\sigma}^2 := \frac{||\hat{\varepsilon}||^2}{n-p}\) (voir ci-dessous à propos de la variance résiduelle).

—- Collapsible section —-

Explications : les résidus

Ecriture non vectorielle : Les résidus sont les différences \(r_i:=y_i-X_i\cdot\hat{\beta}\) entre les valeurs observées de la variable expliquée et le fit associé. Ils sont décorrélés de chacune des variables explicatives.

Ecriture vectorielle : Les résidus, chacun associé à une ligne de \(X\), sont les coordonnées du vecteur suivant :

\[ \hat{\varepsilon} = Y-\hat{Y} = (I-H)\cdot Y = H_{X^\perp}\cdot Y = H_{X^\perp}\cdot\varepsilon. \]

Les résidus appartiennent donc à \(\text{Esp}(X)^\perp\), l’espace dit des résidus. On a que

\[ \begin{aligned} E[\hat{\varepsilon}] & =0 \\ \text{et} \\ Var(\hat{\varepsilon}) &= Var(H_{X^\perp}\cdot\varepsilon) \\ &= H_{X^\perp}\cdot Var(\varepsilon)\cdot H_{X^\perp}^T \\ &= H_{X^\perp}\cdot \sigma^2\cdot I\cdot H_{X^\perp}^T \\ &= \sigma^2\cdot \underbrace{(H_{X^\perp})^2}_{=H_{X^\perp}} \quad\text{(toute matrice de projection est égale à son carré)}\\ &= \sigma^2\cdot H_{X^\perp}^T \end{aligned} \]

On a aussi que

\[ Var(\hat{Y}) = Var(X\cdot\hat{\beta}) = X\cdot Var(\hat{\beta})\cdot X^T = \sigma^2\cdot H \]

—- Collapsible section —-

Explications : la variance résiduelle

Sous l’hypothèse d’homoscédasticité, la variance résiduelle du modèle \(Y=\hat{Y}+\underline{\varepsilon}\) est \(\sigma\cdot I\), i.e. les covariances entre les variables explicatives du modèle sont nulles. Dans ce cas, la variance de la somme de deux variables explicatives est la somme de leurs variances, et la variance globale, appelée aussi variance résiduelle parce que non expliquée par le modèle, est la moyenne des variances de chaque variable explicative. On peut ainsi estimer \(\sigma\) la variance résiduelle par la moyenne des carrés des résidus. Soit \(n\) le nombre de ligne de \(X\) :

\[ \begin{aligned} \text{variance résiduelle } &= \frac{1}{n}\cdot\sum_{i=1}^n\hat{\varepsilon}^2 \\ &=\frac{1}{n}\cdot\sum_{i=1}^n||\hat{\varepsilon}||^2 \\ &= E\left[||\hat{\varepsilon}||^2\right] \\ &= E\left[\text{tr}(\hat{\varepsilon}^T\cdot\hat{\varepsilon})\right] \quad\text{(un scalaire est égal à sa trace)}\\ &= E\left[\text{tr}(\hat{\varepsilon}\cdot\hat{\varepsilon}^T)\right] \\ &= \text{tr}\left(E\left[\hat{\varepsilon}\cdot\hat{\varepsilon}^T\right]\right) \\ &= \text{tr}\left(E\left[\sigma^2\cdot H_{X^\perp}\right]\right) \\ &= \sigma^2\cdot(n-p) \quad\text{pour }p\text{ le nombre de vecteurs colonnes de }X \end{aligned} \]

La dernière égalité provient du fait que la trace d’un projecteur est égale à la dimension du sous-espace sur lequel il projette (voir le théorème de Cochrane). Par conséquent, pour obtenir un estimateur sans biais de \(\sigma^2\), il faut définir : \[ \hat{\sigma}^2 := \frac{||\hat{\varepsilon}||^2}{n-p} \]

Ecriture non vectorielle : L’écart-type est \(\sigma=\sqrt{\sum_{i=1}^n r_i^2/(n-p)}\).

—- Collapsible section —-

Explications : \(R^2\)

Le coefficient de détermination (Bestimmtheitsmass) est noté \(R^2\). Selon Pythagore :

\[ \begin{aligned} ||Y||^2 = ||\hat{Y}||^2+||\hat{\varepsilon}||^2 &= ||\hat{Y}||^2+||\hat{\varepsilon}||^2 \\ &= ||X\cdot\hat{\beta}||^2+||Y-X\cdot\hat{\beta}||^2 \\ \end{aligned} \]

Si la constante fait partie du modèle, donc de \(\text{Esp}(X)\), alors

\[ \begin{aligned} ||Y-\overline{y}\cdot \underline{1}||^2 &= ||\underbrace{Y-\hat{Y}}_{a}+\underbrace{\hat{Y}-\overline{y}\cdot \underline{1}}_{b}||^2 \\ &= ||Y-\hat{Y}||^2+2\cdot\underbrace{(Y-\hat{Y})\cdot(\hat{Y}-\overline{y}\cdot \underline{1})}_{=0}+||\hat{Y}-\overline{y}\cdot \underline{1}||^2 \\ &= ||\hat{Y}-\overline{y}\cdot \underline{1}||^2+||\hat{\varepsilon}||^2 \quad\text{voir schéma}\\ \end{aligned} \]

\[ R^2 := \frac{||\hat{Y}||^2}{||Y||^2} = \cos^2(\theta_0) \quad\text{voir schéma} \]

Si la constante fait partie du modèle, donc de \(\text{Esp}(X)\), on a

\[ \begin{aligned} R^2 &:= \frac{||\hat{Y}-\overline{y}\cdot \underline{1}||^2}{||Y-\overline{y}\cdot \underline{1}||^2} = \cos^2(\theta) \quad\text{, où }\theta\neq\theta_0\text{ , voir schéma} \\ &=1-\frac{||\hat{\varepsilon}||^2}{||Y-\overline{y}\cdot \underline{1}||^2} \end{aligned} \]

Ce coefficient \(R^2\) ne tient cependant aucun compte de la dimension de \(\text{Esp}(X)\), et c’est pourquoi il existe un \(R^2\) ajusté (pour \(p\) le nombre de colonnes de \(X\)) :

\[ R_a^2 := 1-\frac{n-1}{n-p}\cdot\frac{||\hat{\varepsilon}||^2}{||Y-\overline{y}\cdot \underline{1}||^2} \quad\text{si la constante fait partie du modèle} \]

Comment quantifier l’adéquation du modèle avec les données ? Le fit du modèle peut être résumé par \(\hat{\sigma}\) : plus la variance résiduelle est petite, meilleur est le fit. Le fit du modèle peut aussi être résumé par le coefficient de détermination qui est égal à la proportion que prend la variance expliquée par le modèle dans la variance totale. Si on note \(s_y\) l’écart-type des \(y_i\), alors on peut écrire \(r_i=1-\hat{\sigma}^2/s_y^2\). Mais alors deux datasets ayant même droite de régression et même \(\sigma\) mais différant l’une de l’autre par leur étalement sur l’axe des \(y\) (i.e. différant par \(s_y\)) auraient des \(R^2\) différents. Pourtant, leurs fits avec leurs modèles sont évidemment les mêmes puisqu’ils admettent tout deux la même droite de régression. Il est donc difficile d’interpréter \(R^2\).

</div>

summary(fit.3)

## 
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.873 -12.663   2.404  11.356  49.545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.73154    5.87521   4.380 1.49e-05 ***
## mom_hs       5.95012    2.21181   2.690  0.00742 ** 
## mom_iq       0.56391    0.06057   9.309  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.14 on 431 degrees of freedom
## Multiple R-squared:  0.2141, Adjusted R-squared:  0.2105 
## F-statistic: 58.72 on 2 and 431 DF,  p-value: < 2.2e-16

detach(kidiq)

—- Collapsible section —-

Explications : \(t\)-test

On peut vouloir tester si la pente \(\underline{\beta}\) est compatible avec une pente donnée \(\underline{\beta}_0\). La difficulté tient alors au fait que la quantité centrée et normée

\[ \frac{\underline{\hat{\beta}}-\underline{\beta}_0}{\sqrt{Var(\underline{\hat{\beta}})}}\,, \]

dont on pourrait attendre qu’elle suive une loi de Student, ne fait pas sens en l’espèce puisque le dénominateur est une matrice. Par contre on peut tester pour une seule composante du vecteur pente. Soit \(\hat{\beta}_j^0\) le \(j^\text{ème}\) coefficient d’un vecteur \(\hat{\beta}^0\) dont on veut tester si \(\hat{\beta}\) est compatible avec un vecteur de pente donné \(\beta^0\). La statistique du test est alors

\[ T=\frac{\hat{\beta}_j-\hat{\beta}_j^0}{\hat{\sigma}\cdot\sqrt{v_j}}\,, \]

où \(v_j\) est le \(j^\text{ème}\) élément diagonal de la matrice carrée \((X^T\cdot X)^{-1}\). Rappelons à cet effet que \(Var(\underline{\hat{\beta}})=\hat{\sigma}^2\cdot(X^T\cdot X)^{-1}\), et qu’en effet, cette matrice est diagonale en raison du fait que \(Cov(y_i,y_j)=0\) pour \(i\neq j\), par l’hypothèse d’hétéroscédasticité.

En particulier, nous voulons faire ce test pour \(\beta_j^0=0\), ce qui correspond à tester s’il est nécessaire de tenir compte de la \(j^\text{ème}\) variable explicative dans notre modèle de régression, ou s’il vaut mieux ne pas en tenir compte. La statistique devient alors

\[ T=\frac{\hat{\beta}_j}{\hat{\sigma}\cdot\sqrt{v_j}}\,\sim \text{Student à }n-p\text{ degrés de liberté, où }p\text{ est le nombre de vecteurs colonnes de }X. \]

Afin de fixer les idées : Pour une régression simple de la forme \(\beta_0+\beta_1\cdot x+\varepsilon\), on a que \(p=2\).

La commande summary associe à chaque coefficient du vecteur de pente \(\underline{\beta}\) une valeur de cette statistique ainsi qu’une valeur p. Ces deux quantités nous informent quant à l’importance de chacune des variables explicatives dans sa capacité à prédire la variable expliquée. Il faut néanmoins, lorsque \(p\) est grand (autour de 20), préférer la statistique F. Bien sûr la statistique F ne fonctionne que si \(p<<n\). Si \(p>n\), on ne peut même plus ajuster le modèle linéaire aux données et donc on ne peut non plus utiliser cette statistique. Dans ce cas, on a recours à des méthodes de sélection de modèle.

Afin de choisir quelles sont parmi toutes les variables celles qui sont les plus importantes, on peut employer le critère du BIC (Bayesian Information Criterion). On n’oubliera pas, une fois le modèle sélectionné, de calculer la qualité de l’ajustement (le fit) au moyen de, par exemple, le \(R^2\) ajusté.

</div>

—- Collapsible section —-

L’approximation de Laplace

L’approximation de Laplace est nécessaire à comprendre ce qu’est le BIC.

Elle consiste à représenter une distribution postérieure \(p(\theta|D)\) sous forme gaussienne puis à en prendre le début de l’approximation de Taylor. C’est pourquoi nous posons que

\[ p(\theta|D)=\frac{1}{Z}\cdot\exp(-E(\theta))\,, \]

où \(E(\theta)\) est appelée fonction d’énergie, définie comme \(E(\theta):=-\log p(\theta,D)\), et où \(Z:=p(D)\).

Soit \(\theta^*\) le mode de la distribution. L’approximation de Taylor donne

\[ E(\theta)\approx E(\theta^*)+(\theta-\theta^*)\cdot\nabla E(\theta)|_{\theta^*} +\frac{1}{2}(\theta-\theta^*)^T\cdot \left.\frac{\partial^2E(\theta)}{\partial\theta\,\partial\theta^T}\right|_{\theta^*} \cdot(\theta-\theta^*)\,. \]

Or \(\nabla E(\theta)|_{\theta^*}=0\) puisque \(\theta^*\) est le mode de la distribution. Par conséquent, et en notant la hessienne ci-dessus par \(H\) :

\[ \begin{aligned} \hat{p}(\theta,D) &= \exp(-E(\theta)) \\ &= \exp(-E(\theta^*))\cdot\exp\left(-\frac{1}{2}(\theta-\theta^*)^T\cdot H\cdot(\theta-\theta^*)\right) \\ \text{et la distribution postérieure est alors} \\ \hat{p}(\theta|D) &= \frac{1}{Z}\cdot \hat{p}(\theta,D) \quad\propto\mathcal{N}(\theta\,|\,\theta^*,H^{-1}) \\ \text{avec} \\ Z &= p(D) \quad\text{(c'est la vraisemblance)} \\ &= \exp\big(-E(\theta^*)\big)\cdot \left(\sqrt{2\pi}\right)^D\cdot \big(\det(H)\big)^{-\frac{1}{2}} \\ \end{aligned} \]

Cette dernière égalité fournissant une approximation de la vraisemblance. Nous avons en résumé obtenu l’approximation de Laplace de la distribution postérieure et également celle de la vraisemblance. Reécrivons cette dernière approximation (modulo constante multiplicative) :

\[ \log p(D) \approx \underbrace{\log p(D|\theta^*)+\log p(\theta^*)}_{=\log p(\theta^*,D)}-\frac{1}{2}\cdot\log\det(H) \]

On constate que ceci est la \(\log\)-vraisemblance additionnée de deux termes qu’on peut voir comme des termes de pénalité. Si l’on se donne une distribution apriori uniforme, \(p(\theta)\propto 1\), le deuxième terme peut disparaître. Il ne reste alors qu’à substituer \(\theta^*\) par l’estimateur du maximum de vraisemblance \(\hat{\theta}_{MLE}\) pour obtenir

\[ \log p(D) \approx \log p(D\,|\,\hat{\theta}_{MLE})-\frac{1}{2}\cdot\log\det(H) \]

Il faut à ce stade trouver une approximation du dernier terme. C’est pourquoi nous écrivons

\[ \sum_{i_1}^nH_i \quad\text{avec } H_i=\nabla\nabla\log(D_i\,|\,\theta)\,. \]

En approchant les différentes matrices \(H_i\) par une même matrice \(\hat{H}\), on peut écrire

\[ \begin{aligned} \log\det(H) &= \log\det(n\cdot\hat{H}) \\ &= \log \big(n^d\cdot\det(\hat{H})\big) \quad\text{où } d=\dim(\theta)\\ &= d\cdot\log(n)+\log\det(\hat{H})\,. \end{aligned} \]

Or le terme \(\log\det(\hat{H})\) ne dépend pas de \(n\) et est donc nul modulo constante multiplicative. C’est pourquoi nous obtenons l’approximation suivante, dite de Laplace :

\[ \log p(D) \approx \log p(D\,|\,\hat{\theta}_{MLE})-\frac{d}{2}\cdot\log(n)\,. \]

</div>

—- Collapsible section —-

Un critère de sélection de variables : le BIC

Le Bayesian Information Criterion (BIC)

Supposons que nous disposions d’un ensemble de modèles candidats \(\mathcal{M}_m\), \(m=1,2,\ldots,M\) et des paramètres correspondants \(\theta_m\), et que nous souhaitions choisir le meilleur modèle parmi eux. En supposant que nous disposions d’une distribution a priori \(P(\theta_m\,|,\mathcal{M}_m)\) pour les paramètres \(\theta_m\), la probabilité a posteriori d’un modèle donné \(m\) est la suivante

\[ P(\mathcal{M}_m\,|\,Z)=P(\mathcal{M}_m)\cdot P(Z\,|\,\mathcal{M}_m), \]

où \(Z\) représente les données d’apprentissage \(\{(x_i,y_i)\,|\,i=1,\ldots,N\}\). Pour comparer deux modèles \(\mathcal{M}_m\) et \(\mathcal{M}_l\), nous pouvons former les cotes postérieures (the posterior odds)

\[ \frac{P(\mathcal{M}_m)\,|\,Z)}{P(\mathcal{M}_l)\,|\,Z)} =\frac{P(\mathcal{M}_m)}{P(\mathcal{M}_l)}\cdot \frac{P(Z\,|\,\mathcal{M}_m))}{P(Z\,|\,\mathcal{M}_l))}. \]

La quantité \(\frac{P(Z\,|\,\mathcal{M}_m))}{P(Z\,|\,\mathcal{M}_l))}\) est appelée facteur bayésien. Il mesure la contribution des données à la probabilité a posteriori. Si ces probabilités sont supérieures à 1, nous devons choisir le modèle \(m\), le modèle \(l\) dans le cas contraire.

En général, on suppose que l’a priori sur les modèles est uniforme, de sorte que \(P(\mathcal{M}_m)=\text{constant}\). Nous avons besoin d’un moyen d’approcher \(P(Z\,|\,\mathcal{M}_m))\). Nous pouvons utiliser l’approximation de Laplace (voir plus haut) et obtenir

\[ \log P(Z\,|\,\mathcal{M}_m))= \log P(Z\,|\,\hat{\theta}_m,\mathcal{M}_m)) -\frac{d_m}{2}\cdot \log N + \mathcal{O(1)}, \]

où \(\hat{\theta}_m\) est l’estimation ML de \(\theta_m\) et \(d_m\) est le nombre de paramètres libres dans le modèle \(m\).

Si nous choisissons comme fonction de perte \(-2\log P(Z\,|\,\hat{\theta}_m,\mathcal{M}_m))\), nous obtenons le critère d’information bayésien (BIC) sous la forme suivante :

\[ BIC=\frac{N}{\sigma_\epsilon^2}\cdot\left(\overline{err}+\log N\cdot \frac{d}{N}\sigma_\epsilon^2\right), \]

où \(\overline{err}=\frac{1}{\sigma_\epsilon^2}\sum_{i=1}^N(y_i-\hat{f}(x_i))^2\) est l’erreur quadratique moyenne (MSE). Cette expression du \(BIC\) correspond à l’hypothèse que les modèles sont tous gaussiens et que la variance \(\sigma_\epsilon^2\) est connue, alors que la forme générique du \(BIC\) est la suivante

\[ BIC=-2\,\text{loglik}+\log N\cdot d, \] où \(\text{loglik}\) est la \(\log\)-vraisemblance maximisée, c’est-à-dire \(\text{loglik}:=\sum_{i=1}^N\log P_\hat{\theta}(y_i)\) avec \(\hat{\theta}\) l’estimation ML de \(\theta\).

Dans le cadre de la sélection de modèle, choisir le modèle de \(BIC\) minimal correspond à choisir le modèle de plus grande probabilité a posteriori. Et dans le cadre plus particulier de la sélection de variables en régression linéaire, un modèle correspond simplement à une liste de variables sélectionnées, à laquelle on va associer une valeur de \(BIC\).

</div>

attach(kidiq)

fit.2 <- lm (kid_score ~ mom_hs + mom_iq)

plot(mom_iq,kid_score, xlab="Mother IQ score", 
  ylab="Child test score",pch=20, xaxt="n", yaxt="n", type="n")
curve (coef(fit.2)[1] + coef(fit.2)[2] + coef(fit.2)[3]*x, add=TRUE, col="gray")
curve (coef(fit.2)[1] + coef(fit.2)[3]*x, add=TRUE)
points (mom_iq[mom_hs==0], kid_score[mom_hs==0], pch=19)
points (mom_iq[mom_hs==1], kid_score[mom_hs==1], col="gray", pch=19)
axis (1, c(80,100,120,140))
axis (2, c(20,60,100,140))

BIC(fit.2)

## [1] 3768.281

La droite en gris ne considère que les observations associées aux mères qui ont une formation supérieure, tandis que la droite noire ne considère que celles associées à des mères qui n’en ont pas. Les deux droites ont, dans ce modèle, le même coefficient de pente associé au QI de la mère et sont par conséquent parallèles.

—- Collapsible section —-

Inférence dans le modèle gaussien

Jusqu’à ce point, nous n’avons supposé aucun modèle pour décrire les lois des estimateurs. Or si l’on veut pouvoir faire de l’inférence, il faut pouvoir réaliser des tests statistiques qui prennent ces estimateurs comme statistiques de test. Ceci nécessite d’associer une distribution de probabilité aux estimateurs. Précisons d’emblée que certes nous connaissons déjà les variances des estimateurs que nous avons précédemment obtenus mais que que cette connaissance reste insuffisante à conduire un test statistique et proposer des régions de confiance. Pour disposer d’une distribution de probabilité complète pour chaque estimateur, il nous faut faire l’hypothèse d’un modèle. C’est seulement à ce stade que nous allons faire l’hypothèse du modèle gaussien.

Nous donnons ici toutes les hypothèses faites :

\(\text{rang}(X)=p\)
\(E[\varepsilon]=0\) et \(\Sigma_\varepsilon=\varepsilon^2\cdot I\)
La nouvelle hypothèse : \(\varepsilon\sim\mathcal{N}(0,\Sigma_\varepsilon)\)

On remarque que la nouvelle hypothèse implique la deuxième, ainsi que \(\text{Cov}(\varepsilon_i,\varepsilon_j)=\sigma^2\cdot\delta_{ij}\), à savoir que les \(\varepsilon_k\) sont indépendants, mais aussi identiquement distribués de loi \(\mathcal{N}(0,\sigma^2)\).

L’hypothèse faite sur la distribution de \(\underline{\varepsilon}\) nous permet de déduire la distribution des autres estimateurs importants \(\underline{\hat{\beta}}\), \(\hat{\sigma}^2\), \(Var(\underline{\hat{\beta}})\) et \(\hat{Y}\), ce dernier estimateur nous permettant d’obtenir les estimateurs du maximum de vraisemblance pour \(\sigma^2\) et \(\underline{\beta}\).

La vraisemblance de l’échantillon

Nous obtenons ci-après les estimateurs du maximum de vraisemblance pour \(\sigma^2\) et \(\underline{\beta}\). Ce ne sont pas forcément les mêmes que ceux que nous avons déjà obtenus, qui sont les estimateurs des moindres carrés, qu’on a obtenus par des considérations faisant intervenir une matrice de projection.

\[ \begin{aligned} &&L(Y,\underline{\beta},\sigma^2) &= \prod_{i=1}^n f_Y(y_i)& \\ &&&= \left(\frac{1}{2\pi\sigma^2}\right)^{n/2}\cdot\exp\left(-\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n\left(y_i-\sum_{j=1}^p\beta_j\cdot x_i^j\right)^2\right)& \\ &&&= \left(\frac{1}{2\pi\sigma^2}\right)^{n/2}\cdot\exp\left(-\frac{1}{2\sigma^2}\cdot||Y-X\cdot\underline{\beta}||^2\right)& \\ \Rightarrow&&\log L(Y,\underline{\beta},\sigma^2) &= -\frac{n}{2}\cdot\log(\sigma^2)-\frac{n}{2}\cdot\log(2\pi)-\frac{1}{2\sigma^2}\cdot\||Y-X\cdot\underline{\beta}||^2 \\ \end{aligned} \]

Nous obtenons

\[ \begin{aligned} \frac{\partial\log L}{\partial\underline{\beta}} &= \frac{1}{2\sigma^2}\cdot\frac{\partial}{\partial\underline{\beta}}||Y-X\cdot\underline{\beta}||^2 =0 \\ \frac{\partial\log L}{\partial\sigma^2} &= -\frac{n}{2\sigma^2}+\frac{1}{2\pi^4}\cdot||Y-X\cdot\underline{\beta}||^2 =0 \\ \end{aligned} \]

Les solutions de ces équations s’appellent les estimateurs du maximum de vraisemblance, indicés par \(ML\). Ce sont \(\underline{\hat{\beta}}_{ML}=\underline{\hat{\beta}}\) (le même que l’estimateur des moindres carrés) et \(\hat{\sigma}_{ML}^2:=\frac{n-p}{n}\cdot\hat{\sigma}^2\). L’estimateur du maximum de vraisemblance de \(\sigma^2\) est donc biaisé alors que celui des moindres carrés de \(\sigma^2\) ne l’est pas.

Ce sont ces estimateurs, notés génériquement \(\hat{\theta}_{ML}\), que l’on nécessite pour le calcul du BIC, via l’estimation de Laplace (voir ci-dessus).

Lois de ces estimateurs :

\(\underline{\hat{\beta}} \sim \mathcal{N}(\underline{\beta},\sigma^2\cdot(X^T\cdot X)^{-1})\),
\(\hat{\sigma}_{ML}^2 \sim \chi^2_{n-p}\)

Construction d’un intervalle de confiance

Nous souhaitons tester une hypothèse. Pour cela, nous nous donnons un estimateur pertinent. Considérons pour l’exemple l’estimateur d’une moyenne \(\bar{X}=\frac{1}{n}\sum_{i=1}^n x_i\). Nous commençons par centrer et normer cet estimateur : \(T=\frac{\bar{X}-\mu_0}{\hat{\sigma}/\sqrt{n}}\). Nous faisons ceci afin de pouvoir lui associer une loi de probabilité connue. Afin de normer l’estimateur, nous avons besoin de l’écart-type ou d’une estimation \(\hat{\sigma}\) de celui-ci. Nous déterminons ensuite la loi de notre estimateur centré et normé, que nous notons dorénavant \(T\) et appelons notre statistique de test. Dans notre exemple, \(T\) suit une loi de Student à \(n-1\) degrés de liberté.

A partir de l’expression de la statistique \(T\), nous construisons un intervalle de confiance de niveau \(1-\alpha\), où typiquement \(\alpha\) est petit. C’est-à-dire que nous définissons un ensemble de valeurs possibles de \(T\) tel qu’une valeur particulière de \(T\) a seulement une probabilité de \(\alpha\) de ne pas y appartenir. Pour obtenir un tel ensemble, nous le posons égal à \(\left\{T\leq c\right\}\) avec un \(c\) tel que la condition suivante soit satisfaite.

\[ p\left(\left\{T\leq c\right\}\right)=1-\alpha \]

Ceci définit un intervalle de la forme \(]-\infty,b]\) pour une certaine borne \(b\in\mathbb{R}\). Dans notre exemple où l’estimateur à tester est celui d’une moyenne, nous voulons un test bilatéral, i.e. nous prenons comme nouveau test la valeur absolue de l’ancien. Pour le dire autrement, on construit un intervalle à partir de la condition

\[ p\left(\left\{|T|\leq c\right\}\right)=1-\alpha \]

pour obtenir un intervalle de la forme \([a,b]\).

Ci-après nous établissons étape par étape l’intervalle de confiance à partir de cette dernière condition, nous plaçant ainsi dans le cadre d’un test bilatéral.

\[ \begin{aligned} \left\{|\bar{X}-\mu_0|\Big/\frac{\hat{\sigma}}{\sqrt{n}}\leq c\right\} &= \left\{-\bar{X}+\mu_0\leq c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\quad\text{et}\quad \bar{X}-\mu_0\leq c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\} \\ &= \left\{-\bar{X}+\mu_0\leq c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\} \cap \left\{\bar{X}-\mu_0\leq c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\}\\ &= \left\{\bar{X}\geq\mu_0-c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\} \cap \left\{\bar{X}\leq\mu_0+c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\}\\ &= \left\{\mu_0-c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\leq\bar{X}\leq \mu_0+c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right\}\\ &= \left[\mu_0-c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\,,\quad \mu_0+c\cdot\frac{\hat{\sigma}}{\sqrt{n}}\right]\\ \end{aligned} \]

où \(c\) est un quantile donné, par example le 97,5%-tile de la loi de la statistique de test (ici, une loi de Student à \(n-1\) degrés de liberté) si le test est bilatéral. Notons que cet intervalle de confiance est bien de niveau \(1-5\%=95\%\) et non de niveau \(1-5\%/2=97,5\%\). En effet,

\[ p\left(\left\{|\bar{X}-\mu_0|\Big/\frac{\hat{\sigma}}{\sqrt{n}}\leq c\right\}\right)=95\%\,. \]

</div>

—- Collapsible section —-

Erreur de prédiction

Prédiction

Considérons une nouvelle observation, non déjà contenue parmi les \(n\) lignes de \(X\) la matrice des données. C’est un vecteur de dimension \(p=\)nombre de vecteurs-colonnes de \(X\), vecteur que nous notons \(x_{n+1}=(x_{n+1}^1,\ldots,x_{n+1}^p)^T\). Nous voulons prédire la valeur associée \(y_{n+1}\).

Nous savons que \(y_{n+1}=x_{n+1}\cdot\underline{\beta}+\varepsilon_{n+1}\),

où \(E[\varepsilon_{n+1}]=0\), \(Var(\varepsilon_{n+1})=\sigma^2\) et \(Cov(\varepsilon_{n+1},\varepsilon_i)=0\quad\forall i=1,\ldots,n\).

La prédiction est donc \(\hat{y}_{n+1}=x_{n+1}\cdot\underline{\hat{\beta}}\).

Dans cette expression, \(x_{n+1}\) est une grandeur déterministe donnée et donc la seule source de variabilité provient de \(\underline{\hat{\beta}}\). Par conséquent l’erreur \(\varepsilon_{n+1}\) sur la prédiction est simplement \(x_{n+1}\) multiplié par l’erreur sur \(\underline{\hat{\beta}}\).

Or l’estimateur du vecteur de pente

\[ \underline{\hat{\beta}}=\frac{||\hat{\varepsilon}||^2}{n-p}\cdot(X^T\cdot X)^{-1} \]

est le produit de \(||\hat{\varepsilon}||^2\) avec un facteur déterministe, si bien que les lois des estimateurs \(\underline{\hat{\beta}}\) et de \(||\hat{\varepsilon}||^2\) sont toutes deux de la même famille, en l’occurence de loi de Student, une fois ceux-ci centrés et normés. Or nous connaissons la variance de \(\underline{\hat{\beta}}\), ce qui nous permet justement de normer l’estimateur \(\underline{\hat{\beta}}\).
Bon… Dit ainsi, c’est faux parce que \(\underline{\hat{\beta}}\) est une grandeur vectorielle et sa variance, une matrice tandis que \(||\hat{\varepsilon}||^2\) est une quantité scalaire. Tout de même, le raisonnement tient si l’on considère, à la place de ce vecteur, ses composantes scalaires \(\hat{\beta}_i\). (Voir à ce propos les explications sur le \(t\)-test ci-dessus).

L’erreur de prédiction \(y_{n+1}-\hat{y}_{n+1}\)

Nous calculons la variance de l’erreur de prédiction \(y_{n+1}-\hat{y}_{n+1}=y_{n+1}-x_{n+1}\cdot\underline{\hat{\beta}}\).

Nous notons que l’erreur de prédiction n’est pas l’erreur sur la prédiction ! La seconde est la partie non déterministe du modèle décrivant la quantité prédite \(\hat{y}_{n+1}\) tandis que la première est la différence entre les deux quantités aléatoires que sont l’observation \(y_{n+1}\) et la quantité prédite. C’est à présent celle-ci qui nous intéresse.

\[ \begin{eqnarray*} Var(y_{n+1}-\hat{y}_{n+1}) &=& Var(x_{n+1}\cdot\underline{\beta}+\varepsilon_{n+1} -x_{n+1}\cdot\underline{\hat{\beta}}) \\ &=& \sigma^2+x_{n+1}^T\cdot Var(\underline{\hat{\beta}})\cdot x_{n+1} \\ &=& \sigma^2\cdot\left(1+x_{n+1}^T\cdot\left(X^T\cdot X\right)^{-1}\cdot x_{n+1}\right) \end{eqnarray*} \]

C’est l’incertitude sur l’observation \(y_{n+1}\), i.e. \(\sigma^2\), additionnée de l’incertitude sur \(\hat{y}_{n+1}\). C’est bien sûr également \(Var(\hat{\varepsilon}_{n+1})\). Cette expression permet de définir un intervalle de confiance et ainsi d’associer graphiquement à la droite de régression une bande d’erreur de prédiction.

Intervalle de confiance bilatéral de niveau (\(1-\alpha\)), pour l’erreur de prédiction d’une nouvelle observation \(y_{n+1}\) : \[ \left[ x_{n+1}^T\cdot\hat{\beta}\,\pm\, t_{n-p}(1-\alpha/2)\cdot \hat{\sigma}\cdot \sqrt{1+x_{n+1}^T\cdot(X^TX)^{-1}\cdot x_{n+1}}\right], \]

où \(t_{n-p}(1-\alpha/2)\) est le fractile de niveau \(1-\alpha/2\) de la loi de Student à \(n-p\) degrés de liberté.

La façon détaillée d’obtenir un tel intervalle est expliquée ci-dessus, sous “Inférence dans le modèle gaussien”, à partir d’un exemple.

Quant à comprendre pourquoi cette statistique suit une loi de Student, nous faisons les considérations suivantes sur la statistique de test \(T\) dont découle notre intervalle de confiance :

\[ T=\frac{y_{n+1}-\hat{y}_{n+1}}{\hat{\sigma}\cdot\sqrt{1+x_{n+1}^T\cdot(X^TX)^{-1}\cdot x_{n+1}}} \]

\(\underline{\hat{\beta}}\) suit une loi normale et \(x_{n+1}\) est fixe donc \(\hat{y}_{n+1}\) suit une loi normale. La valeur aléatoire à prédire \(y_{n+1}\) suit une loi normale \(\mathcal{N}(x_{n+1}\cdot\underline{\hat{\beta}},\,\sigma^2)\) et est indépendante des \(y_1,\ldots,y_n\) par hypothèse. Nous avons donc que \(y_{n+1}\) est indépendante de \(\hat{y}_{n+1}\) puisque \(\hat{y}_{n+1}\) est une fonction linéaire des \(y_1,\ldots,y_n\). L’erreur de prédiction, i.e. le numérateur de \(T\), suit par conséquent une loi normale. Et dans la mesure où l’écart-type au dénominateur est estimé et non pas connu, la fraction suit bien une loi de Student.

Remarque : Le propos, lorsque l’on qualifie \(y_{n+1}\) de “nouvelle” observation, est que cette observation n’a absolument pas servi pour le calcul des estimateurs ni de \(\underline{\hat{\beta}}\) ni de \(\hat{\sigma}^2\). En ce sens elle est nouvelle, et c’est ce que son indice \(n+1\) suggère, étant bien clair que la matrice \(X\) ne contient que \(n\) lignes, et que donc les estimateurs ont été ajustés avec \(n\) observations, dont la “nouvelle” ne fait pas partie.

</div>

Ajoutons pour chacune des droites de régression les intervalles de confiance pour \(\hat{\beta}\) et pour l’erreur de prédiction. On affiche tout d’abord les paramètres du modèle :

attach(kidiq)
fit.2 <- lm (kid_score ~ mom_hs + mom_iq)
fit.2[["coefficients"]]

## (Intercept)      mom_hs      mom_iq 
##   25.731538    5.950117    0.563906

coef(fit.2)[1]

## (Intercept) 
##    25.73154

A l’aide de cette information sur le modèle, nous adaptons la formule des droites de régression en modifiant les deux commandes curve(), considérant qu’une fois la variable mom_hs vaut 0 et qu’une fois elle vaut 1 :

plot(mom_iq,kid_score, xlab="Mother IQ score", 
  ylab="Child test score",pch=20, xaxt="n", yaxt="n", type="n")
curve (coef(fit.2)[1] + coef(fit.2)[2] + coef(fit.2)[3]*x, add=TRUE, col="red", lwd = 3)
curve (coef(fit.2)[1] + coef(fit.2)[3]*x, add=TRUE,lwd = 3)
points (mom_iq[mom_hs==0], kid_score[mom_hs==0], pch=19)
points (mom_iq[mom_hs==1], kid_score[mom_hs==1], col="grey", pch=19)

grille <- seq(min(kidiq[,"mom_iq"]),max(kidiq[,"mom_iq"]),length=100)
grille <- data.frame(grille)

# Ajout des deux intervalles de confiance à la droite associée à mom_iq=0
vector_of_zeros <- rep(0, 100)
grille_0 <- cbind(grille,vector_of_zeros)
grille_0 <- grille_0[,c("vector_of_zeros","grille")] #swapping the columns
colnames(grille_0) <- c("mom_hs","mom_iq")
ICdte_0 <- predict(fit.2,new=grille_0,interval="confidence",level=0.95)
ICprev_0 <- predict(fit.2,new=grille_0,interval="pred",level=0.95)
matlines(grille_0$mom_iq,cbind(ICdte_0,ICprev_0[,-1]),lty=c(1,2,2),col=1)

# Ajout des deux intervalles de confiance à la droite associée à mom_iq=1
vector_of_ones <- rep(1, 100)
grille_1 <- cbind(grille,vector_of_ones)
grille_1 <- grille_1[,c("vector_of_ones","grille")] #swapping the columns
colnames(grille_1) <- c("mom_hs","mom_iq")
ICdte_1 <- predict(fit.2,new=grille_1,interval="confidence",level=0.95)
ICprev_1 <- predict(fit.2,new=grille_1,interval="pred",level=0.95)
matlines(grille_1$mom_iq,cbind(ICdte_1,ICprev_1[,-1]),lty=c(1,2,2),col=2)

axis (1, c(80,100,120,140))
axis (2, c(20,60,100,140))

Modèles avec interactions

En ajoutant une interaction entre les variables explicatives, le modèle fera deux droites non plus forcément parallèles pour chacune des deux catégories de niveau de formation des mères.

fit <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
display (fit)

## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
##               coef.est coef.se
## (Intercept)   -11.48    13.76 
## mom_hs         51.27    15.34 
## mom_iq          0.97     0.15 
## mom_hs:mom_iq  -0.48     0.16 
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.23

plot(mom_iq,kid_score, xlab="Mother IQ score", 
  ylab="Child test score",pch=20, xaxt="n", yaxt="n", type="n")
curve (coef(fit)[1] + coef(fit)[2] + (coef(fit)[3] + coef(fit)[4])*x, add=TRUE, col="gray")
curve (coef(fit)[1] + coef(fit)[3]*x, add=TRUE)
points (mom_iq[mom_hs==0], kid_score[mom_hs==0], pch=20)
points (mom_iq[mom_hs==1], kid_score[mom_hs==1], col="gray", pch=20)
axis (1, c(80,100,120,140))
axis (2, c(20,60,100,140))

BIC(fit)

## [1] 3765.451

Ajoutons pour chacune des droites de régression les intervalles de confiance pour \(\hat{\beta}\) et pour l’erreur de prédiction. On affiche tout d’abord les paramètres du modèle, qui inclut ici un terme d’interaction :

attach(kidiq)
fit.2 <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
fit.2[["coefficients"]]

##   (Intercept)        mom_hs        mom_iq mom_hs:mom_iq 
##   -11.4820211    51.2682234     0.9688892    -0.4842747

coef(fit.2)[1]

## (Intercept) 
##   -11.48202

plot(mom_iq,kid_score, xlab="Mother IQ score", 
  ylab="Child test score",pch=20, xaxt="n", yaxt="n", type="n")
curve (coef(fit.2)[1] + coef(fit.2)[2] + coef(fit.2)[3]*x + coef(fit.2)[4]*x, add=TRUE, col="red", lwd = 3)
curve (coef(fit.2)[1] + coef(fit.2)[3]*x, add=TRUE,lwd = 3)
points (mom_iq[mom_hs==0], kid_score[mom_hs==0], pch=19)
points (mom_iq[mom_hs==1], kid_score[mom_hs==1], col="grey", pch=19)

grille <- seq(min(kidiq[,"mom_iq"]),max(kidiq[,"mom_iq"]),length=100)
grille <- data.frame(grille)

# Ajout des deux intervalles de confiance à la droite associée à mom_iq=0
vector_of_zeros <- rep(0, 100)
grille_0 <- cbind(grille,vector_of_zeros)
grille_0 <- grille_0[,c("vector_of_zeros","grille")] #swapping the columns
colnames(grille_0) <- c("mom_hs","mom_iq")
ICdte_0 <- predict(fit.2,new=grille_0,interval="confidence",level=0.95)
ICprev_0 <- predict(fit.2,new=grille_0,interval="pred",level=0.95)
matlines(grille_0$mom_iq,cbind(ICdte_0,ICprev_0[,-1]),lty=c(1,2,2),col=1)

# Ajout des deux intervalles de confiance à la droite associée à mom_iq=1
vector_of_ones <- rep(1, 100)
grille_1 <- cbind(grille,vector_of_ones)
grille_1 <- grille_1[,c("vector_of_ones","grille")] #swapping the columns
colnames(grille_1) <- c("mom_hs","mom_iq")
ICdte_1 <- predict(fit.2,new=grille_1,interval="confidence",level=0.95)
ICprev_1 <- predict(fit.2,new=grille_1,interval="pred",level=0.95)
matlines(grille_1$mom_iq,cbind(ICdte_1,ICprev_1[,-1]),lty=c(1,2,2),col=2)

axis (1, c(80,100,120,140))
axis (2, c(20,60,100,140))

Interprétation des coefficients

fit <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
display (fit)

## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
##               coef.est coef.se
## (Intercept)   -11.48    13.76 
## mom_hs         51.27    15.34 
## mom_iq          0.97     0.15 
## mom_hs:mom_iq  -0.48     0.16 
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.23

Il y a une certaine difficulté à interpréter les coefficients dès lors que le modèle admet un terme d’interaction. Par exemple, le coefficient de 51,27 est la différence de QI que le modèle prédit entre des enfants dont la mère n’a pas de formation supérieure (mom_hs=0) et ceux dont la mère a une formation supérieure (mom_hs=1), parmi tous les enfants dont la mère a un QI de zéro. Et puisque de telles mères n’existent pas, le coefficient de 51,27 est sans signification réelle. De même, le coefficient de 0,97 est l’augmentation en QI que prédit le modèle pour chaque augmentation d’une unité de QI de la mère parmi les enfants dont la mère n’a pas de formation supérieure. Dans ce cas-ci, ce coefficient signifie quelque chose.

On peut toutefois simplifier l’interprétabilité des coefficients en centrant chacune des variables explicatives :

attach(kidiq)
c.mom_hs <- mom_hs - mean(mom_hs)
c.mom_iq <- mom_iq - mean(mom_iq)
fit <- lm (kid_score ~ c.mom_hs + c.mom_iq + c.mom_hs:c.mom_iq)
display (fit)

## lm(formula = kid_score ~ c.mom_hs + c.mom_iq + c.mom_hs:c.mom_iq)
##                   coef.est coef.se
## (Intercept)       87.64     0.91  
## c.mom_hs           2.84     2.43  
## c.mom_iq           0.59     0.06  
## c.mom_hs:c.mom_iq -0.48     0.16  
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.23

Cette fois, le coefficient de 2,84 s’interprète comme la différence prédite entre les moyennes des deux groupes (enfants dont la mère a ou n’a pas de formation supérieure) parmi les enfants dont la mère a un QI moyen mean(mom_iq) (=100). Quant au coefficient de 0,59, il prédit une augmentation de 0,59 unités de QI par unité de QI supplémentaire de la mère parmi les enfants dont la mère est mean(mom_hs) (=0,79). Tout cela fait davantage sens, n’est-ce pas ? On peut bien évidemment ne centrer que la variable mom_iq et garder inchangée mom_hs, pour une interprétation vraiment plus sensée.

Le centrage des variables explicatives aide ainsi à interpréter les coefficients de pente. Cependant le coefficient de la variable catégorielle mom_hs est bien plus grand que celui de la variable métrique mom_iq. Ceci parce que l’on compare le changement complet entre l’une et l’autre de deux catégories de formation avec le changement d’une seule unité de QI de la mère. Pour changer cela, il est possible de diviser chacune des variables explicatives par deux écarts-type.

z.mom_hs <- (mom_hs-mean(mom_hs))/(2*sd(mom_hs))
z.mom_iq <- (mom_iq-mean(mom_iq))/(2*sd(mom_iq))
fit <- lm (kid_score ~ z.mom_hs + z.mom_iq + z.mom_hs:c.mom_iq)
display (fit)

## lm(formula = kid_score ~ z.mom_hs + z.mom_iq + z.mom_hs:c.mom_iq)
##                   coef.est coef.se
## (Intercept)       87.64     0.91  
## z.mom_hs           2.33     1.99  
## z.mom_iq          17.65     1.82  
## z.mom_hs:c.mom_iq -0.40     0.13  
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.23

On choisit de diviser par non pas par un mais par deux écarts-type afin de maintenir la cohérence quand on considère une variable binaire (0/1). On voit cela en considérant une simple variable binaire \(x\) valant soit 0 soit 1 avec probabilité égale \(p=0.5\). L’écart-type est alors \(\sqrt{0.5\cdot 0.5}=0.5\) et ainsi la variable standardisée \((x-\mu_x)/(s\cdot\sigma_x)\) prend les valeurs \(\pm 0.5\), et son coefficient reflète la comparaison entre \(x=0\) et \(x=1\).

Par contraste, si l’on avait divisé par un seul écart-type, la variable standardisée prendrait les valeurs \(\pm 1\) et son coefficient refléterait la moitié de la différence entre les deux valeurs possibles des \(x\). Le raisonnement tient aussi pour des variables binaires avec \(p\neq 0.5\) en raison du fait que \(\sqrt{p(1-p)}\approx0.5\) lorsque \(p\) n’est pas trop différent de 0,5.

La régression polynomiale

attach(kidiq)
fit <- lm(kid_score ~ poly(mom_iq, 3))
fit[["coefficients"]]

##      (Intercept) poly(mom_iq, 3)1 poly(mom_iq, 3)2 poly(mom_iq, 3)3 
##         86.79724        190.39123        -61.15731         29.81028

plot(mom_iq,kid_score, xlab="Mother IQ score", 
  ylab="Child test score",pch=20, xaxt="n", yaxt="n", type="n")
curve (coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 + coef(fit)[4]*x^3, add=TRUE, col="red", lwd = 3)
points (mom_iq[mom_hs==0], kid_score[mom_hs==0], pch=19)
points (mom_iq[mom_hs==1], kid_score[mom_hs==1], col="grey", pch=19)

grille <- seq(min(kidiq[,"mom_iq"]),max(kidiq[,"mom_iq"]),length=100)
grille <- data.frame(grille)

# Ajout des deux intervalles de confiance à la courbe polynomiale
vector_of_zeros <- rep(0, 100)
grille_0 <- cbind(grille,vector_of_zeros)
grille_0 <- grille_0[,c("vector_of_zeros","grille")] #swapping the columns
colnames(grille_0) <- c("mom_hs","mom_iq")
ICdte_0 <- predict(fit,new=grille_0,interval="confidence",level=0.95)
ICprev_0 <- predict(fit,new=grille_0,interval="pred",level=0.95)
matlines(grille_0$mom_iq,cbind(ICdte_0,ICprev_0[,-1]),lty=c(1,2,2),col=2)

axis (1, c(80,100,120,140))
axis (2, c(20,60,100,140))

Les résidus

Intéressons-nous au plot des résidus pour le modèle le plus simple, à un seul prédicteur. Ce plot ne dévoile aucune tendance qui pourrait suggérer une analyse en série temporelle. On le voit et le test du portemanteau conduit ci-après le confirme.

## Fit the model
fit.2 <- lm (kid_score ~ mom_iq)
resid <- fit.2$residuals
sd.resid <- sd(resid)

 # Figure 3.12
plot (mom_iq, resid, xlab="Mother IQ score", ylab="Residuals", pch=20)
abline (sd.resid,0,lty=2)
abline(0,0)
abline (-sd.resid,0,lty=2)

Nous testons la blancheur des résidus en calculant la statistique de Ljung-Box.

require(caschrono)
ret=c(1,2,3,4,5,10,15,20)
# les différentes valeurs de retard choisies pour le test du portemanteau
Box.test.2(resid,nlag=ret,type="Ljung-Box",decim=2)

##      Retard p-value
## [1,]      1       0
## [2,]      2       0
## [3,]      3       0
## [4,]      4       0
## [5,]      5       0
## [6,]     10       0
## [7,]     15       0
## [8,]     20       0

La régression logistique

La régression logistique est la façon habituelle de modéliser les variables à expliquer binaires (binary outcomes). Nous commençons par modéliser la préférence politique (vote républicain : 1, vote démocrate : 0) lors des élections présidentielles de 1992 en fonction du salaire des votants, salaire discrétisé en valeurs entières allant de 1 (pauvre) à 5 (riche).

Nous rappelons à toutes fins utiles que

\[ p(y_i=1)=\text{logit}^{-1}(X_i\cdot\underline{\beta})\,, \]

où \(X_i\cdot\underline{\beta}\) est appelé le prédicteur linéaire, \(X_i\) une ligne de la matrice \(X\), et \(\text{logit}^{-1}(x)=\frac{\exp(x)}{1+\exp(x)}\).

nes <- read_dta('data/nes/nes5200_processed_voters_realideo.dta')
attach(nes)
yr <- 1992
  ok <- year==yr & presvote<3
  vote <- presvote[ok] - 1
  income <- income[ok]
head(nes)

## # A tibble: 6 × 62
##    year resid weight1 weight2 weight3 age   gender  race    educ1   urban region
##   <dbl> <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl+l> <dbl+l> <dbl+l> <dbl> <dbl+>
## 1  1948  1001       1       1       1 NA    1 [1. … 1 [1. … 1 [1. … NA    NA    
## 2  1948  1002       1       1       1 NA    2 [2. … 1 [1. … 2 [2. … NA    NA    
## 3  1948  1003       1       1       1 NA    2 [2. … 1 [1. … 2 [2. … NA    NA    
## 4  1948  1004       1       1       1 NA    2 [2. … 1 [1. … 3 [3. … NA    NA    
## 5  1948  1005       1       1       1 NA    1 [1. … 1 [1. … 3 [3. … NA    NA    
## 6  1948  1006       1       1       1 NA    2 [2. … 1 [1. … 2 [2. … NA    NA    
## # ℹ 51 more variables: income <dbl+lbl>, occup1 <dbl+lbl>, union <dbl+lbl>,
## #   religion <dbl+lbl>, educ2 <dbl+lbl>, educ3 <dbl+lbl>,
## #   martial_status <dbl+lbl>, occup2 <dbl+lbl>, icpsr_cty <dbl+lbl>,
## #   fips_cty <dbl+lbl>, partyid7 <dbl+lbl>, partyid3 <dbl+lbl>,
## #   partyid3_b <dbl+lbl>, str_partyid <dbl+lbl>, father_party <dbl+lbl>,
## #   mother_party <dbl+lbl>, dlikes <dbl+lbl>, rlikes <dbl+lbl>,
## #   dem_therm <dbl+lbl>, rep_therm <dbl+lbl>, regis <dbl+lbl>, …

fit.1 <- glm (vote ~ income, family=binomial(link="logit"))
display(fit.1)

## glm(formula = vote ~ income, family = binomial(link = "logit"))
##             coef.est coef.se
## (Intercept) -1.30     0.18  
## income       0.30     0.06  
## ---
##   n = 1222, k = 2
##   residual deviance = 1623.5, null deviance = 1655.0 (difference = 31.5)

Interpretation des coefficients

On considère la relation suivante :

\[ \text{logit}\Big(p(y_i=1\,|\,X_i)\Big)=X_i\cdot\underline{\beta}\,, \]

qui nous fait comprendre qu’un changement d’une unité d’une des composantes \(x_i\) du vecteur \(X_i\) fait changer le prédicteur de \(\beta_i\) unités, mais le changement de \(\beta_i\) unités de l’expression \(\text{logit}(\ldots)\) n’est pas constant mais dépend de \(X_i\). En l’espèce, une interprétation du changement d’une unité de l’un des prédicteurs \(x_i\) n’est pas possible. Par contre on peut linéariser la fonction \(\text{logit}^{-1}\) autour de la moyenne. La dérivée de \(\text{logit}^{-1}(\alpha+\beta x)\) par rapport à \(x\) (mis ici pour le prédicteur \(x_i\)) donne

\[ \frac{\partial}{\partial x}\text{logit}^{-1}(\alpha+\beta x) = \frac{\beta\cdot\exp(\alpha+\beta x)}{(1+\exp(\alpha+\beta x))^2}\,, \]

et par conséquent la pente au point où le prédicteur \(x\) prend la valeur centrale \(\bar{x}\) vaut

# Compute the mean over the non-NA values
mean_income <- mean(income, na.rm = TRUE)
slope_at_mean <- fit.1$coef[2]*exp(fit.1$coef[1]+fit.1$coef[2]*mean_income)/
  (1+exp(fit.1$coef[1]+fit.1$coef[2]*mean_income))^2
slope_at_mean

##     income 
## 0.07326868

fit.1$coef[2]/4

##     income 
## 0.07582499

On constate qu’il suffit grosso-modo de diviser \(\beta\) par 4 pour obtenir la pente en \(\bar{x}\). Et puisque la pente en ce point est maximale, on peut interpréter la quantité \(\beta_i/4\) comme une majoration de l’augmentation de la probabilité \(p(y_i=1)\) pour une augmentation d’une unité du prédicteur \(x_i\).

On peut aussi interpréter les coefficients en termes de “odds ratios”. En effet, on a que

\[ \begin{aligned} \log\left(\frac{p(y=1\,|\,x)}{p(y=0\,|\,x)}\right) &= \alpha+\beta x \\ \text{On peut essayer de s'en convaincre :} \\ \frac{p(y=1\,|\,x)}{p(y=0\,|\,x)} &= \exp(\alpha+\beta x) \\ \frac{\text{logit}^{-1}(X_i\cdot\underline{\beta})}{1-\text{logit}^{-1}(X_i\cdot\underline{\beta})} &= \exp(\alpha+\beta x) \\ \frac{\frac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)}}{1+\frac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)}} &= \exp(\alpha+\beta x) \\ \frac{1}{1+\frac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)}} \cdot \frac{1}{1+\exp(\alpha+\beta x)} &= 1 \\ \frac{1}{1+\exp(\alpha+\beta x)+\exp(\alpha+\beta x)} &= 1+\exp(\alpha+\beta x) \\ \frac{1}{1+2*\exp(\alpha+\beta x)} &= 1+\exp(\alpha+\beta x) \\ \end{aligned} \]

Une augmentation d’une unité du prédicteur \(x\) correspond à une multiplication de l’odd ratio par le facteur \(\exp(\beta)\). Mais ce concept d’odd ratio reste obscure, et c’est pourquoi nous nous en tiendrons à la précédente interprétation en terme de majoration de l’augmentation de probabilité.

Graphiquement :

# Graph figure 5.1 (a)
 curve (invlogit(fit.1$coef[1] + fit.1$coef[2]*x), 1, 5, ylim=c(-.01,1.01),
         xlim=c(-2,8), xaxt="n", xaxs="i", mgp=c(2,.5,0),
         ylab="Pr (Republican vote)", xlab="Income", lwd=4)
  curve (invlogit(fit.1$coef[1] + fit.1$coef[2]*x), -2, 8, lwd=.5, add=T)
  axis (1, 1:5, mgp=c(2,.5,0))
  mtext ("(poor)", 1, 1.5, at=1, adj=.5)
  mtext ("(rich)", 1, 1.5, at=5, adj=.5)
  points (jitter (income, .5), jitter (vote, .08), pch=20, cex=.1)

Pour voir un exemple d’analyse de données menée de bout en bout et faisant usage de la régression logistique, voir pp.86-104 du livre de Andrew Gelman: Data Analysis Using Regression and Multilevel/Hierarchical Models, 2006, Cambridge University Press.